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Abstract 


The article deals with constructing a Toeplitz-like preconditioner for linear 
systems arising from finite difference discretization of the spatial fractional 
diffusion equations. The coefficient matrices of these linear systems have an 
S+L structure, where S is a symmetric positive definite (SPD) matrix and 
L satisfies rank(L) < 2. We introduce an approximation for the SPD part 
S, which is called Ps, and then we show that the preconditioner P = Ps +L 
has the Toeplitz-like structure and its displacement rank is 6. The analysis 
shows that the eigenvalues of the corresponding preconditioned matrix are 
clustered around 1. Numerical experiments exhibit that the Toeplitz-like 
preconditioner can significantly improve the convergence properties of the 
applied iteration method. 
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1 Introduction 


In this article, we aim to propose a Toeplitz-like preconditioner for solving 
the linear system resulting from a finite difference approximation of an initial- 
boundary value problem of the spatial fractional diffusion equations (FDEs) 
that were introduced in [4]. In this article, the authors show that applying 
the finite difference method to the FDEs leads to the following linear system: 


Ax = b, (1) 
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where the coefficient matrix of the linear system (1) is a dense Toeplitz-like 
matrix. As we see in [4], the coefficient matrix of (1) can be rewritten as the 
sum of a symmetric positive definite matrix S and a low rank matrix L, that 
is, A= 5 +L (for more details, see [4, 2]). The matrix S has the following 
structure: 


S =nIn+ GD,GT + G7 DG, (2) 


where 7 > 0, G is a Toeplitz matrix, and D, and Dz, are diagonal matrices 
with positive diagonal elements (more details are available in Section 2). In 
[4], the authors proposed a preconditioner by replacing G with its Strang’s 
circulant approximation and two matrices D, and Dz by their scalar (identity 
matrix) approximation. They showed the superlinear convergence of their 
preconditioner, when D, and Dz are scalar identity matrices. In [2], we 
introduced a preconditioner based on the approximation of S' defined in (2), 
and we proved that the proposed preconditioner is strong even though D, 
and Dz are not a multiple of the identity matrix. In the following, we give 
some definitions about Toeplitz-like matrices. 

Toeplitz-like matrices are defined by means of a displacement operator. 
Given an n X n matrix A, we consider the down-shift matrix of order n, 
ie (e2 €3.°°' En On), and the displacement operator V defined by 


V(A) =A-Z,AZP. (3) 


We define the displacement rank of matrix A as rank(V(A)). The matrix 
A is said to have the Toeplitz-like structure if its displacement rank is low 
compared with its order n, that is, rank(V(A)) =m <n. Ifthe displacement 
rank of A is m, then we can rewrite V(A) as V(A) = CD’, where C,D € 
R”"*™. Two matrices C and D are called the generators of the matrix A. 

Throughout this article, we use the following notations: Capital letters, 
boldface lowercase letters, and regular lowercase letters denote matrices, vec- 
tors, and scalars, respectively. Moreover, J, denotes the identity matrix of 
order n, while J, denotes the nxn exchange matrix J, = antidiag(1,1,...,1). 
We denote by e; the jth column of identity matrix, and 0 denotes the zero 
vector of proper dimensions. 

The organization of this article is as follows. In Section 2, we present 
the new preconditioner for the discretized linear systems arising from FDEs. 
Some computational remarks are given in Section 3. Numerical experiments 
are provided in Section 4 to prove the performance of our preconditioner. 
Finally, in Section 5, some concluding remarks are provided. 


2 Preconditioning technique 


The finite difference discretization of FDEs was developed in [4]. By intro- 
ducing 
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the authors showed that the coefficient matrix of resulting linear system by 
finite difference discretization of FDEs possesses A = S+L € R"*” structure, 
where 


S=nn+ GD, G@? + GT pD_G™, (5) 
L= cag! + eJa(Jg)’, (6) 


in which 7 > 0, D4 = diag(din,di.2,...,din) with d-; > 0, 1,c2 > 0, 
and 


Nes eee 

BY GOP see 

Gi) = m ee € RX”, (7) 
eg a 


For simplicity, in the rest of the the superscript, (a) will be ignored. We see 
that the matrix S' is a symmetric positive definite and that D is a matrix with 
low rank, that is, rank(L) < 2. Denote by dimin and dimax the smallest 
and the largest elements of Di, respectively. Letting d+ = Seman tO min in 
[2], we used the following preconditioner : 


Pg = nln + d4GG? + d_G'G. (8) 


This preconditioner is an approximation of S. We can replace S in A by Ps, 
obtaining so-called TL preconditioner 


P=Ps+L (9) 
for the matrix A. 


Theorem 1. If P is defined as (9), then P~'A = M+N, where M is similar 
to P5*S and rank(N) < 2. 


Proof. We can extend S = span{g, Jg} to an orthonormal basis of R”. Let 
V= (vi Vo-°° Vn) be the matrix representation of such orthonormal basis, 


T 
=. and v2 = Eine where h = — GEE. Suppose that V 
2 2 


such that v, = 


98 Akhoundi 


is partitioned as V = (Vi V2), where V; = (v1 V2) and V2 = (v3 Vae0° Va) 
From (6), we see that LV2 = 0. We have 


(10) 


T T p-l 


Ve P5'LV, 0 
If we define kK = Ig + Vi Ps *LVi and 
_ Ih-K "VF — O*V Tp 
Q =V (sean we 0 Vv ’ 
then 
(I+ P51L)1 =V [I+ V7 PsILV] Vt 


V 
-1 

7 Vii Ps LV, 0 T 

=v [+ (rein) 


K- 0 \ Ur 
=v (seo -) . 
I+ 


We note that QV2 = 0. By this assumption, we can prove our assertion in 
the following 


P1A=(Pg+L)(S+L)=(1+ Ps'L)'P5'S(I+ SL) 
= (I+ P5'L)-1Ps'S(1+ Pg L)\(I+ Py*L)-'(1+S8""L). (11) 


In (11), if we define M = (I + Pg'L)-1P5*S(I + P5'L), then 


P'A=M(I+Q)I+S7'L) 
=M+N, 


where N = M(S-'L+Q+QS7'!L). We see that M is similar to Pes, It 
can be easily verified that NV2 = 0. Hence rank(N) < 2. 


Theorem 2. [2] If \ € o(P;'S), then |1 — \| < k < 1, where k = aH <1 


: a n a di jmaxtdyjmin d—,maxt+d—,min 
with 7 = ae and « = max{ 3 ; 5 \. 
Based on Theorems | and 2, the eigenvalues of the preconditioned matrix 
P-'A are clustered around one, and thus Krylov subspace methods with the 


proposed preconditioner coverage very quickly. 


Remark 1. The main differences between our proposed preconditioner P 
and the proposed preconditioner in [2] Ps are as follows: 


1- Pgs is an approximation of S in (5), while P consists of two parts. The 
first part Ps is an approximation of the matrix S and the second one 
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is the matrix L (the rest of the matrix A). Hence, we expect that P is 
more accurate than Ps. 


2- In the next section, in order to show the computational efficiency of pre- 
conditioner P, we prove that P has the Toeplitz-like structure, and its 
generators will be computed precisely. By using them, we can use fast 
algorithms to compute P~'. Hence, in comparison with the circulant 
approximation of Ps in [2], the P~' is computed straightforward. 


In Section 4, numerical experiments show that the efficiency of the proposed 
preconditioner P. 


3 Some computational results 


In this section, we show that the matrix P defined in (9) has the Toeplitz- 
like structure, and we construct its generators. To this end, we define the 
following temporary vectors: 


= T : T 
& = (90 91 -** Jn-1) f= (91 g2°*+Gn-1) - (12) 
The following relations can be easily deduced from (4) and (12) 


—ao — (a1 a ag) 


1 —a1 — (a2 — a1) 
= = = Zr n&n> 13 
a-—g T@—a) at ane (13) 
—aAn-1 — (dn _ Gn—1) 
gs = ge + Zg, (14) 
Z's = 8— Gn€n = (5) ; (15) 


We use the following lemmas in our subsequent discussion. 


Lemma 1. Let G be the n x n matrix defined in (7). Then the following 
relations hold: 


(i) V(GG*) = gg’, 
(ii) V(GTG) = e18"G + G’gef — Bg’ geref — JZ7 gg" ZJ, 
(iii) V(Ps) = nee? +d, ge? +d_(e,g7 G+G? gel —g’ ge,e) —JZ gg" ZJ). 


Proof. The matrix G can be partitioned as 


Gx ( ae (16) 
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where G; = G(2 : n,2 : n) is a lower triangular Toeplitz matrix, so the 
displacement structure of GG" can be viewed as 
ViGG")=GG = 266 7 
~\t Gi 0 Gt G,0/\0 0 )° 


Now, by a straightforward computation, part (7) of the lemma can be verified 


For part (iz), first we know that 
(18) 


eG = (ge tC). 
Therefore, 
V(G?G) = G"G-— za" azt (19) 
_ (9 t™\ (9 0) (0 0\ (0 G (20) 
0OGT/\tG GE IEP AOE 
ols +r 
eg tG ) 
= ae aes 21 
(f8 —JttT J e? 


0 


It gives that 
e’g 0\ (gg 0\) (0 0 
0 0 Jt td 


Spe. “5p 
Toy (8 8 Gi 
vera) = (83 0 )+ (Es § 0 
=e.8'G4+G" gel — ge gele, — JZ" ge" ZI. (22) 
We note that V(I) = ee? , so 
V(Ps) = nee? +. d,.V(GG") +d. =W(G"@) 
= neie; +d,8g" +d_(eig’G+G" ge; —g’ geile, — JZ" gB"ZJ), 


which completes the proof of (ii) 


Lemma 2. If we define the auxiliary vectors hy = d_(Jg) — c2Z(Jg) and 
he = cig — d4(Zg), then the following relations hold: 


(i) ¢ 2 [(Ja)(Je)" — Z(Ja)(Jg)? Z7|— a JE" ee" 2d = (c2—d_)(Ja)(Jg)"+ 
Z(Jah? + e1(d_ay + d_gn)(Jg)* + d_(9n(Jg) — greet, 
(dy. — c1)(Za)(Zg)” + ahy + 


(ii) cr [(a)(g)” — Z(a)(g)7Z7] + dy eg? = 
e1(aod4 + d4go)(Zg)" + (d+90(Zg) + dygxe1)el. 


Proof. We see that 
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co [(Ja)(Jg)” — Z(Ja)(Jg)" Z" | — d_JZ7 gg" ZI 
= ¢2 [(Ja)(Jg)” — Z(Ja)(Jg)" 27] — d_ [J(g — gnen)(@—gnen)* J] By (15) 
= [e9(Ja) — d_(Jg)] (Jg)* — e2Z(Ja)(Jg)" Z* 


+ d_gne1(Jg)” + d_gn(Jg)et — d_gheie, 
= J [cga— d_(a— Z™a—anen)| (Jg)" —e2Z(Ja)\(Jg)"Z" ~—_ By (13) 
+ d_gnei(Jg)" +d_gn(Je)et — d_gneret 


= (co — d_)(Ja)(J g)* + Z(Ja)hy + e1(d_ An + d_ Ga ey” 
+ (d_gn(Jg) _ d_g2ei)eq. 


In the similar way, we can prove part (i) as follows: 


1 [(a)(g - Za ene + d,gg7 

= c1 [(a)( (a)(g)"Z"] + ds(goe1 + Zg)(goe1 + Ze)" By (14) 

=Z(dig— ary (Zg)? + cag’ +d [9oe1e, + goei(Zg)? + 9028! | 
=Z [(d. Gt ee AGS d4An€n| (Zg)? + cag” By (13) 


+d, [goe1et + goe1(Zg)" + 90Zge} | 
(dy. — c1)(Za)(Zg)" + ahd + e1(aody + dyg0)(Zg)* 
+ (d4.go(Zg) + dy.ggei)e, 


Theorem 3 indicates that P in (9) has the Toeplitz-like structure, and we 
can use a superfast solver to compute P~!. 


Theorem 3. The matrix P defined in (8) is a Toeplitz-like matrix and its 
displacement operator is 


V(P) =(c2 — d_)(Ja)(Jg)? + (dy — c1)(Za)(Zg)" + Z(Ja)hy 
+ah3 +e,h] + hie?, (23) 
where h, and hg are defined in Lemma 2, and 


hg = (d_an + d_gn)(Jg) + d+(ao + go)Zg + ne: + d_G*g, 


ha = d_gn(Jg) — gne1 + d4.90(Zg) + dygge1 + d_-G"g — d_g" ger. 


Proof. We have P = Ps + L. Hence 


V(P) =V(Ps) + V(L) = V(Ps) + e1((a)(g)” — Z(a)(Zg)*) 
+ e(Ja(Jg") — ZJa(ZJg)*). (24) 


By part (ii) of Lemma 1, we can rewrite (24) as follows: 
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V(P) = neiet + d,eg" + d_(eig’G + G’get — g’geiet — JZ" gg" ZJ) 
+e1((a)(g)" — Z(a)(Zg)") + co(Ja(Jg") — ZJa(ZJg)") 
= c2 [(Ja)(Jg)" — Z(Ja)(Jg)" Z| — d_JZ" gg" ZI 

+c; [(a)(g)? — Z(a)(g)? 27] + d, 8a” 
+neie, + d_(eig'G+G" get — g’ Bele, ) 
Parts (z) and (iz) of Lemma 2 and the above equations imply that 


V(P) =(c2 — d_)(Ja)(Jg)* + Z(Ja)hy + e1(d_an + d_gn)(Je)” 
+ (Lgn(Je) - ger)e? + (dy — 01)(Za)(Ze)? + ah? 
+e; (aod + dy90)(Zg)* + (dsgo(Zg) + dygoer)e, 
+neret + d_(eig’G + G’ get — g gee, ) 
=(c2 — d_)(Ja)(Jg)" + Z(Ja)ht + (dy — c1)(Za)(Zg)" + ah 
+ ey ((d_an 1 d_gn)(Jg) EW d4(ao + 90)Zg + Nei + d_G™g)" 
+ (d_9n(Jg) — gner + d490(Zg) + dy.goe1 + d_G? g — d_g* Bex) ey 
=(cz — d_)(Ja)(Jg)” + (dy — c1)(Za)(Zg)” 
+ Z(Ja)yhi +ah5 + e.h? + hie? . 


We showed that our preconditioner P has the Toeplitz-like structure and 
its displacement rank of P (rank(V(P))) is 6. So we can use fast and sta- 
ble Levinson like algorithms [5] to compute P~'r in O(n?) operations. For 
Krylov subspace iteration methods such as preconditioned GMRES [10], high 
quality preconditioning plays a crucial role in accelerating the convergence 
speed of the Krylov subspace iteration methods. Instead of O(n”) operations 
to apply our preconditioner in each iteration, Circulant preconditioners re- 
quire only O(n log(n)) operations in each step. Significant reduction of total 
iterations (and total CPU time) by using our preconditioner in numerical 
experiments show the efficiency of our preconditioner. 


4 Numerical experiments 


The GMRES(20) (restarts every 20 inner iterations) with the proposed pre- 
conditioner is applied to solve the linear system (1). We choose the right-hand 
side such that «* = (1,2...,n)" is the exact solution, that is, b = Ax*. The 
stopping criterion in the numerical experiments is ||r,||2/||b||2 << le—7, where 
rz is the residual vector of the linear system after k iterations and b is the 
right-hand side. For all experiments, the initial guess is chosen as the zero 
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Table 1: Numerical results for Example 1 
a n I Po Ps TL 
Iter CPU Iter CPU Iter CPU Iter CPU 
0.7] 27 | 81 23.5113] 6 0.6189 | 4 0.4012 | 1 0.3001 
qu T - 7 2.1718 5 1.6312 1 0.9017 
ae || - 7 7.9142 5 5.1305 1 2.7321 
Q18 T - 8 9.8320 5 6.3483 1 5.9821 
ore |) : 11 14.8639) 7 9.6407 1 8.9231 
0.9 | 2% A . 13 1.4921 8 1.2134 1 0.5743 
gu T - 21 9.5901 13 5.8513 1 1.2401 
gle ai - 29 46.2163 | 17 26.2245 1 3.5421 
2 : t : t - 2 7.3142 
qi4 T = t 7 tT = 2 10.6893 


vector. All the numerical experiments are run in MATLAB on a desktop with 
the configuration: Intel(R) Core(TM) CPU Q9450 2.66 GHz and 8.00 GB 
RAM. In the following tables, “J” represents the GMRES method without 
preconditioning technique. We use Po to denote the circulant preconditioner 
[4], Ps to denote preconditioner in [2], and TL to present the proposed pre- 
conditioner defined in (9). Also the “Iter” denotes the number of iterations, 
“CPU” denotes the total CPU time in seconds for solving the problem, and n 
denotes the size of our linear system. In the tables, { indicates no convergence 
within 100 iterations. 


Example 1. In this example, we consider 


1 2. 
D. =di a PP ait ol 
+ diag {0, (——) =) , ’ i (25) 
ae nm—2.,_, ,2-38\1_¢ 
D_ = diag{1,("—) =z? pear Os (26) 


and cy = cp = 2. 


Example 2. In this example, D+ are random diagonal matrices with diag- 
onal entries taken independently at random from [0,1] and c; = cg = 2. 


The numerical results of Examples 1 and 2 are shown in Tables 1 and 2, 
respectively. Tables 1—2 show clearly the fast convergence of the TL precon- 
ditioner. We see that T'L preconditioner is more effective in both CPU time 
and the number of iterations than Ps and Po preconditioners. 

In what follows, we compare the eigenvalues of A and P~!A for Example 
1, the results are shown in Figures 1 and 2. As we see the eigenvalues of 
preconditioned matrix P~!A are clustered around 1. 
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Table 2: Numerical results for Example 2 
Qa n I Po Pg TL 
Iter CPU | Iter CPU Iter CPU Iter CPU 
0.7/2 | 94 0.5931] 6 0.6907 4 0.4455 1 0.2517 
que} + - 7 1.7012 5 1.3259 1 0.8993 
ge} + : 7 8.2302 5 6.4216 1 2.3263 
7 aa i | a 9 11.4134] 6 8.1983 1 7.3129 
a |, - 11 143921 7 11.0214) 2 10.9386 
0.9 | 2° | + - 30 2.3102 | 17 ~~ 1.7492 1 0.4352 
que} + : 41 14.0561} 25 10.9318 | 1 1.3107 
oO .|) 75 47.5169 | 44 20.9312} 1 4.3563 
oe | : t - t : 2 9.8311 
qi} + 3 + z t x 2 14.0968 
1.5 1.5 T 
1F 1 1 7 
0.57 al 0.5 4 


0G@0OGD0 O O OC OCD 


Figure 1: Eigenvalues of P~1A (left) and A (right) for Example 1, with a = .9 and 


ae 


5 Concluding remarks 


In this article, the preconditioned GMRES method with Toeplitz-like struc- 
ture was employed to solve the discretized linear systems arising from FDEs. 
The displacement structure of the Toeplitz-like preconditioner was computed 
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9 fe r r 1.5 r 1 


0.5 1 1.5 0 50 100 


Figure 2: Eigenvalues of P~1A (left) and A (right) for Example 1, with a = .9 and 
n= 21, 


precisely. The efficiency of the proposed preconditioner was proved even 
though the diffusion coefficients are not constants. Numerical experiments 
have demonstrated the efficiency of the proposed preconditioner. 

Nevertheless, it is interesting that to propose the possible Toeplitz-like 
preconditioner for two-dimensional FDEs, and then examine the efficiency of 
the proposed preconditioner. 
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